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We performed numerical simulations on a two-dimensional inelastic hard disk system under 
gravity with a heat bath to study the dynamics of granular fluidization. Upon increasing the 
temperature of the heat bath, we found that three phases, namely, the condensed phase, locally 
fluidized phase, and granular turbulent phase, can be distinguished using the maximum packing 
fraction and the excitation ratio, or the ratio of the kinetic energy to the potential energy. It is 
shown that the system behavior in each phase is very different from that of an ordinary vibrating 
bed. 
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The dynamics of fluidized granular systems has at- 
tracted much attention in physics, .communities as a 
nonequilibrium statistical systemBH'B Due to the fact 
that an element particle in granular systems is already 
macroscopic, there are two major differences in their dy- 
namics from that of an ordinary molecular system. First, 
thermal fluctuation does not play any role because rel- 
evant energy scales for the kinetic and potential energy 
are much larger than the thermal energy. Second, the 
dynamics is dissipative because the degrees of freedom 
we are dealing with are coupled with microscopic pro- 
cesses that are not treated explicitly. Because of these 
properties, granular systems require an energy source 
in order to be in a steady state. It is also expected 
that the external gravitational force plays an impor- 
tant role in their dynamics. In most experimental situa- 
tions where the systems are excited Jpitiaj/ibrating plate, 
the effect of gravity is very large .& l is 1313' B It has been 
found that various patterns and localized oscillations ap- 
pear as the, iutmsity and frequency of the vibration are 
changed S3 It has been also demonstrated that the 
velocity distribution of the particles does not have the 
Maxwell-Boltzmann distribution—reflecting the fact that 
the system is not in equilibrium. LSItirlla' Systems devoid 
of gravity are also being examined mainly, .by; c.ompiJtcr 
simulations and theoretical analysesOOESIij'EJ'LilL^ 
A one-dimensional system with a heat bath at one end 
has been investigated and it has been found that the 
spatial distribution of particles is singular, therefore, the 
hydrodynamical description of the behavior does not ap- 
pear to be possibleEj) Systems with uniform excitation, 
in which all the particles _axe agitated by the Langevin 
force, are also studiedEJEiP because these system are 
statistically homogeneous and have a thermodynamic 
limit, even though such a situation may not be easily 
realized experimentally. Large density fluctuation and 
non-Gaussian velocity distribution have been observed 
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in such systemsciP when the effect of particle collision is 
larger than that of the Langevin force. Since the grav- 
itational force should have a large effect on clustering, 
it would be of interest to observe how the system be- 
havior changes as the gravity sets in, which is one of 
the motivating factors behind this work. An experimen- 
tal study for analyzing the relative effects of gravity to 
those of excitation has been reported.E3 In the exper- 
iment, stainless steel balls were placed on a plate that 
was held almost horizontally, and were excited by the 
vibrating wall of the lower side. By tilting the system, 
the effective gravity was changed and it was shown that 
even a small force of gravity makes the cluster migrate 
downward. In this letter, we examine the effect of grav- 
ity on granular fluidization by numerical simulations. It 
is demonstrated that a system undergoes a few phase 
changes as the ratio of the external driving to the grav- 
ity is changed. 

The system considered here consists of hard disks of 
mass m and diameter d in a two-dimensional space un- 
der uniform gravity. We employ the periodic boundary 
condition in the horizontal direction and the system is 
considered to be high enough to prevent the particles 
from hitting the ceiling too frequently. The particles col- 
lide with each other with a restitution constant r. All 
the disks are identical, namely the system is monodis- 
persed. For simplicity, we neglect the rotational degree 
of freedom. Because of the hard core interaction, col- 
lision is instantaneous and only binary collisions occur. 
When two disks, i and j, with respective velocities 
and Uj collide, the velocities after the collision, and 
u'- , are given by 

u- = Ui - -(1 + r)[n- ( Ui - Uj)]n (1) 

Uj = u, + i(l + r)[n ■ (m - u^n, (2) 

where n is the unit vector parallel to the relative position 
of the two colliding particles in contact. Between the 
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colliding events, the particles undergo free fall motion 
with the gravitational acceleration g = (0, —g) following 
parabolic trajectories. The system is driven by a heat 
bath with temperature T w at the bottom of the system; 
a disk hitting the bottom bounces back with velocity v = 
(v x , v y ) chosen randomly by the probability distributions 
<t>x(v x ) and <t>y{v v )\ 



2nkfiT, 



_ e -mvi/2k B T w ^ <Vx< 
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-v v e 
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(0 <v y < oo), (4) 



where k B is the Boltzmann constant. e-P We assume that 
the collision with the ceiling is elastic with an unchanged 
horizontal velocity. Note that the wall temperature T w 
here is just a parameter to characterize the external driv- 
ing and is not related to the thermal fluctuation. We call 
this system "inelastic hard disk system with a heat bath 
under weak gravity" (IHSHG) to emphasize the impor- 
tance of the competition between excitation of the heat 
bath and gravity. The system is completely characterized 
with only four dimensionless parameters; the restitution 
coefficient r, the driving intensity A = ksT w /mgd, the 
system width in the unit of disk diameter N w = L/d, 
and the number of layers Nh = N/N W) where L is the 
system width and N is the total number of disks. The 
present system-is analogous to the ordinary granular vi- 
brated bedjclLT&tj'E^F but simpler because it does not 
have an external time scale, or a period of external driv- 
ing vibration. We employ the event-driven molecular dy- 
namics method to simulate the IHSHG model; wp have 
developed a simple and efficient algorithm,Ej'L3' Elf 1 and 
achieved the fastest simulation speed in the world (about 
460 million collisions per CPU hour for the 4000 disk 
system on a VT-Alpha-600), which allows us to simulate 
the system over a wide range of parameters. Most of the 
simulations presented in this letter were performed on 
the system with r = 0.9, N w = 200, and N h = 20 (i.e. 
N — 4000). The driving intensity A was varied from 110 
to 770 to study the change in the system behavior. To 
ensure that the system had reached a steady state, we 
relaxed the system until the total energy did not drift. 
Thereafter, we started the simulations with the length of 
10,000 collisions per particle to obtain the data. Inelastic 
collapseEJ did not occur in the present system during the 
simulation time for the parameter region studied here. 

In Fig. 1, a typical snapshot of the simulation is shown 
for A = 182 with (r,N w ,N h ) = (0.9,100,20). The av- 
erage packing fraction profile as a function of height y 
is also shown. It can be seen that the density is small 
near the bottom because of the excitation by the heat 
bath, and the packing fraction reaches the maximum 
value A max at a certain height. 

In order to characterize the steady states, we measure 
the maximum packing fraction A max as a function of the 
driving intensity A (o in Fig. 2). There are two cusps 
around A = 200 and 380, suggesting phase transitions 
with changing A. These transitions should be related to 
the excitation structure of the state, therefore we define 
the excitation ratio [i by [i = K./U, namely, the ratio of 
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Fig. 1. A typical snapshot of the simulation and the average 
packing fraction profile as a function of height y. The parameters 
are (r, A, N w ,N h ) = (0.9, 182, 100, 20). The black disks are the 
disks in the layer where the packing fraction reaches a maximum 
value. 



the total kinetic energy K. = m/2 v 1 to the potential 
energy U = mgJ^iUi- The cusps appear at the same 
points for which confirms the transitions (x in Fig. 2). 
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Fig. 2. The maximum packing fraction A ma x and the excitation 
ratio fi vs. the driving intensity A. The system parameters are 
(r, N w ,N h ) = (0.9, 200, 20). The regions of A for the condensed 
phase, the locally fluidized phase, and the granular turbulent 
phase are shown. 



In order to conceptualize the underlying physical 
mechanism of the system behavior in each phase sep- 
arated by the above transitions, three snapshots for typ- 
ical values of A for each phase are shown in Fig. 3. For 
A = 182 (al-a3 in Fig. 3), there is a closed packing layer 
and the state is not very dynamic, because excitation is 
low and the potential energy is dominant. We call the 
phase in A < 200, the condensed phase. It is apparent, 
however, that collective motion appears on the surface 
of the layer. In the second phase, A = 250 (bl-b3 in 
Fig. 3), the closed packing layer is locally broken by ex- 
citation of the heat bath. The high-speed particles are 
blown upward from the holes in the layer. In Fig. 3, we 
see only one hole in the system, but for a larger system, 
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Fig. 3. The series of snapshots for the three typical values of A, 
al-a2 (A = 182), bl-b2 (A = 250), cl-c2 (A = 667). The system 
parameters are (r, N w , N h ) = (0.9, 200, 20). 



there are certain cases where we can observe more than 
one hole. The holes migrate and occasionally become 
more active; sometimes they almost close temporarily, 
but the structure is fairly stable. We call this phase 
the locally fluidized phase. For the case of A = 667 
in the third phase (cl-c3 in Fig. 3), the layer is com- 
pletely destroyed. The average density is quite low, but 
it is very different from the ordinary molecular gas phase. 
The density fluctuation is very large and this fluctuation 
causes turbulent motion driven by the gravity. At some 
time, the whole system gets excited with some relatively 
smaller density fluctuations, but the very next moment, 
a large proportion of the particles travels downward and 
forms a layer like structure. This structure, however, is 
destroyed immediately. We call this phase the granular 
turbulent phase. 

These inhomogeneous behaviors should result in a non- 
Maxwell-Boltzmann distribution of velocity. Figure 4(a) 
shows the horizontal velocity distribution functions at 
the horizontal layer around the height of the maximum 
packing fraction for each phase in the log-linear scale. 
The distributions for A=250 and 667 deviate from the 
Gaussian and are more or less in exponential form in 
the tail region, ft can be seen that the distribution for 
A = 250 in the locally fluidized phase consists of two 
parts; the central part that originates from the closed 
packing layer, and the wide tail that comes from the 
fluidized holes. The central part is very close to that for 
A = 182. 

In order to quantify the deviation, we calculated the 
flatness parameter / =< v^. > / < v\ > 2 (Fig. 5). Over 
the whole region, the value of / is different from 3, which 
is the value for the Gaussian distribution, but it is re- 
markable that / becomes very large, as large as 20, in 
the locally fluidized phase. This unexpected large value 
of /, however, can be understood naturally as follows. 
Assume the PDF cf>(v) has two components, the nar- 
row Gaussian distribution 4>g{v) with the weight 1 — p 
and the broader stretched exponential distribution <ps (v) 
with the weight p; 

<t>(v)=p4> s (v) + (l-p)4> G (v), (0<p<l), (5) 
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Fig. 4. The horizontal velocity probability distribution functions 
for the three typical values of A. The other parameters are 
(r, N w ,N h ) = (0.9, 200, 20). fa) The data are plotted in the unit 
of yjk^fjm. (b) Equation (g) is fitted to the data for A = 250 
with the parameters uq = 0.048, xo/crc = 1.7, /3 = 0.85, and 
p = 0.25. Note that the data are scaled by the standard devia- 
tion a. 

where cf> s (x) = (P/\2x T(l/ f3)])e^ x / Xo ^ and <j) G (x) = 
[l/(V2^cr G )]e" :E2/2 ' TG . Then, the flatness of this distri- 
bution is given by 

[r(5//3)/r(l/ j 3)]^p + 3(l-p) _ ~ _ %o_ (6) 

• ([r(3/ / 3)/r(i//3)]^ P + (i- P )) 2 ' °g 

because the second and the fourth moment of 4>s{x) are 
given by T(3/0)/T(l/0) ■ x\ and T(5/ (3)/Th/p) ■ 4, re- 
spectively. If we fit the parameters in eq. (gj) to the data 
for A = 250, we obtain a G = 0.048, i = 1.7, (1 = 0.85, 
and p = 0.25, which yields / ~ 19.0 (Fig. 4(b)). The 
enhancement of the flatness by the superposition of the 
broader distribution can even be drastic for a small value 
of p as can be seen if eq. (j^) is plotted as a function of p. 
This observation indicates that / can be used as a sen- 
sitive index to detect the appearance of a small weight 
of the broad component in the distribution. The sharp 
rise of / in Fig. 5 around A = 200 is clear evidence of 
the appearance of the fluidized holes. 

There are some interesting issues regarding the dy- 
namics in each phase. 

(1) What is the mechanism of the surface wave-like 
motion in the condensed phase? There must be no sur- 
face tension in the granular system, therefore the restor- 
ing force should originate from the balance between the 



1 



Masaharu Isobe and Hiizu Nakanishi 



20 



10 



Condensed ; 


Locally 
Fluldlzed 


Granular Turbulent 




% 




k 


m 

• 


i L 






* 






Exponential 
Gaussian 



'o 200 400 600 800 



&+ 

Fig. 5. The flatness / vs. the driving intensity A. The other 
parameters are (r,N w ,N h ) = (0.9,200,20). The values for the 
Gaussian distribution (/ = 3) and the exponential distribution 
(/ = 6) are indicated in the figure. The statistical error bars 
with twice the standard deviation are also plotted in the several 
typical data for each phase. 



gravity and the excitation, but it is not clear if it could 
be described as an ordinary gravitational wave in fluid. 

(2) The localized excitation in the locally fluidized 
phase should resemble a^drcle in a 3-d system and re- 
minds us of an oscillonjij) but they are different; the 
external vibration frequency is essential for the oscillon 
dynamics, but the localized excitation here does not have 
such a characteristic frequency. In bl-b3 of Fig. 3, we 
can observe the merging process of two excitations, but 
their mode of interaction is not clear. An interesting is- 
sue here is that if the transition from the fluidized phase 
to the condensed phase and/or the transition to the gran- 
ular turbulent phase can be understood in terms of the 
local excitation, does the locally fluidized phase trans- 
form to the condensed phase when the distance between 
the excitations diverge? Does the locally fluidized phase 
transform to the granular turbulent phase at the point 
where they condense? 

(3) Turbulent motion does not exist in the 1-d hard 
rod system.EiP It should be very different from ordinary 
fluid turbulence. In the finite system with finite height, 
turbulence may disappear when the driving A is large 
enough, but in the system with infinite height, it appears 
that turbulent motion persists however large the driving 
A is. 

Before concluding, let us discuss the relationship be- 
tween the present model and the ordinary vibrated bed, 
where the system is driven by a vibrating plate at a given 
frequency. In such a case, the control parameter is usu- 
ally taken as the ratio of the accelerations T = Aoj 2 /g, 
where A and u are the amplitude and the angular fre- 
quency of the vibration, respectively. In most experi- 
mental situations, the external frequency dominates the 
system time scale and the collision interval is directly 
related to 1/cj. In the present case, the control param- 
eter is the ratio of the energies A = ksT w /mgd and no 



external time scale is imposed. This situation may cor- 
respond to that in the vibrated bed when l/u> is much 
shorter than any other of the relevant time scales in the 
system. 

In summary, we performed numerical simulations for 
a two-dimensional inelastic hard sphere system with a 
heat bath under gravity. Upon increasing the heat bath 
temperature, the system exhibited three distinct phases, 
namely, the condensed phase, the locally fluidized phase, 
and the granular turbulent phase. Their dynamics are 
very different from those of not only an ordinary molec- 
ular system but also of a conventional vibration bed sys- 
tem. 
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